Turing instability in quantum activator–inhibitor systems

Turing instability is a fundamental mechanism of nonequilibrium self-organization. However, despite the universality of its essential mechanism, Turing instability has thus far been investigated mostly in classical systems. In this study, we show that Turing instability can occur in a quantum dissipative system and analyze its quantum features such as entanglement and the effect of measurement. We propose a degenerate parametric oscillator with nonlinear damping in quantum optics as a quantum activator–inhibitor unit and demonstrate that a system of two such units can undergo Turing instability when diffusively coupled with each other. The Turing instability induces nonuniformity and entanglement between the two units and gives rise to a pair of nonuniform states that are mixed due to quantum noise. Further performing continuous measurement on the coupled system reveals the nonuniformity caused by the Turing instability. Our results extend the universality of the Turing mechanism to the quantum realm and may provide a novel perspective on the possibility of quantum nonequilibrium self-organization and its application in quantum technologies.

www.nature.com/scientificreports/ state was discussed [67][68][69] , which can be regarded as a quantum manifestation of the Turing-type bifurcation originally analyzed in a classical system 70 . Though this bifurcation is interesting, it is not exactly the Turing instability in the original sense because the considered system is not of the activator-inhibitor type and does not possess a homogeneous stationary state when the coupling is absent 70 . Additionally, the relation between the Turing bifurcation and quantum features, such as quantum entanglement and quantum measurement, has not been studied in these papers [67][68][69] .
In this study, we analyze Turing instability in the original sense of Turing 23 and Gierer and Meinhardt 24 in quantum dissipative systems in the simplest setting, i.e., in a pair of symmetrically coupled units, by providing a minimal model of quantum activator-inhibitor systems. We show that a degenerate parametric oscillator with nonlinear damping can behave as a quantum activator-inhibitor unit and that diffusive coupling between two such units can induce Turing instability and lead to nonuniformity and entanglement between the two units, which gives rise to a pair of nonuniform states that are symmetrically mixed due to quantum noise. We further demonstrate that performing continuous measurement on the coupled system breaks this symmetry and reveals the true asymmetry caused by the Turing instability. A schematic diagram is shown in Fig. 1.

Quantum activator-inhibitor system
Quantum activator-inhibitor unit. We first show that a single-mode, degenerate parametric oscillator with nonlinear damping in quantum optics 71 can be considered a quantum activator-inhibitor unit in the sense that the deterministic trajectory of the system in the classical limit obeys conventional activator-inhibitor dynamics.
We denote by ω 0 the resonance frequency of the cavity and by ω p the frequency of the pump beam of squeezing. In the rotating coordinate frame of frequency ω p /2 , the evolution of the density operator ρ representing the system state obeys the quantum master equation (QME) 71 where [A, B] = AB − BA is the commutator of two operators A and B, a is the annihilation operator that subtracts a photon from the system, a † is the creation operator that adds a photon to the system ( † denotes the Hermitian conjugate), � = ω 0 − ω p /2 is the detuning of the resonance frequency of the system from the half frequency of the pump beam, ηe iθ ( η ≥ 0 ) is the squeezing parameter representing the effective amplitude of the pump beam, D[L]ρ = LρL † − (ρL † L − L † Lρ)/2 is the Lindblad form representing the coupling of the system with the reservoirs through the operator L ( L = a or L = a 2 ), and γ 1 (> 0) and γ 2 (> 0) are the decay rates for linear and nonlinear damping, i.e., the single-photon and two-photon loss, respectively, due to coupling of the system with the respective reservoirs. The reduced Planck constant is set as = 1.
We employ the phase-space method 72,73 and use the Wigner distribution W(x, p) as the quasiprobability distribution to represent the density operator ρ , where x and p denote the position and momentum in the phase space, respectively. Using this approach, we can transform the QME to the evolution equation for W(x, p) on the phase space, which generally has derivative terms higher than the second order. When γ 2 is small, we can neglect the higher order derivative terms, and the evolution equation for W(x, p) corresponding to QME (1) can be approximated by a semiclassical Fokker-Planck equation (FPE) or the corresponding stochastic differential equation (SDE). The deterministic trajectory in the classical limit of QME (1), which neglects the effect of small quantum noise and is given by the deterministic part of the SDE, is found to obey the following two-dimensional system: See "Methods" for the detailed derivation of the equations and characterization of the quantum regime. By appropriately choosing the parameters, classical system (2) obeys activator-inhibitor dynamics (see "Methods"). We set the parameters such that the position x and momentum p play the roles of the activator and inhibitor variables, respectively, namely, x autocatalytically enhances its own production while p suppresses the growth of x. It is noted that the system without nonlinear damping can also behave as a quantum activator-inhibitor unit, but nonlinear damping is necessary to prevent the system state from diverging to infinity after destabilization at the origin. Figure 2a shows the deterministic vector field of Eq. (2), where the two curves represent nullclines of x and p (on which ẋ = 0 or ṗ = 0 ) and their intersection at (x, p) = (0, 0) corresponds to a stable fixed point. Figure 2b shows a scatter plot of a single trajectory of the semiclassical SDE obtained by direct numerical simulations (DNSs) in the steady state (see "Methods"), and Fig. 2c shows the stationary Wigner distribution obtained from QME (1). The semiclassical trajectory and the Wigner distribution are distributed around the classical fixed point at the origin due to quantum noise.
Diffusively coupled quantum activator-inhibitor units. In the classical Turing instability, the uniform stationary state of spatially distributed activator-inhibitor systems is destabilized when diffusion of the activator and inhibitor species with appropriate diffusivity is introduced, leading to the formation of nonuniform states 23 . In the simplest setting, this counterintuitive Turing instability can already be observed in a system consisting of two diffusively coupled activator-inhibitor units with identical properties: a uniform stationary state of the system, in which the two units take the same states, becomes destabilized when the diffusivities are appropriately chosen, resulting in the formation of a nonuniform stationary state, in which the two units settle into different states from each other. As a quantum model that undergoes Turing instability, we diffusively couple two identical quantum activator-inhibitor units (denoted 1 and 2), each of which obeys Eq. (1). The coupled system of the two units is described by a two-mode density operator ρ , which obeys the QME where a j and a † j are the annihilation and creation operators for the jth quantum activator-inhibitor unit ( j = 1, 2 ), respectively. The parameters �, ηe iθ , γ 1 and γ 2 are common to both units. In this equation, the first line represents the two single-mode units given by Eq. (1), and the newly introduced terms in the second line represent the coupling between the two units. The first coupling term can be represented as a sum of squeezing terms, i.e., , ρ , which can be interpreted as single-mode and two-mode squeezing Hamiltonians, respectively. The second term with D c represents dissipative coupling, namely, a coupling arising from dissipative processes 12,14 . It is noted that Eq. (3) is symmetric with respect to the exchange of the units 1 and 2.
By employing the phase-space method for two-mode systems, the deterministic dynamics in the classical limit of QME (3) can be derived as (see "Methods") . www.nature.com/scientificreports/ where x j and p j represent the position and momentum of the jth unit in the phase space of the two-mode Wigner distribution W(x 1 , p 1 , x 2 , p 2 ) 73 . We see that two classical activator-inhibitor units, each of which is described by Eq.
(2), are diffusively coupled through the position x (activator) and momentum p (inhibitor) by the last term in each equation. These terms arise from the single-and two-mode squeezing Hamiltonians whose intensities are characterized by D h and from the dissipative coupling whose intensity is characterized by D c in Eq. (3). The diffusion constants of x and p in Eq. (4) are given by D x = (D c + D h )/2 and D p = (D c − D h )/2 , respectively. It should be noted that the first term characterized by D h represents a Hamiltonian coupling and non-dissipative, but it acts as a dissipative coupling in the deterministic dynamics in the classical limit in Eq. (4). The classical coupled system described by Eq. (4) can undergo Turing instability when the conditions of local self-enhancement and long-range inhibition are satisfied (see "Methods"). Therefore, the quantum activator-inhibitor system, Eq. (3), is also expected to exhibit Turing instability when the parameter values are appropriately chosen. Our aim in this study is to clarify whether Turing instability can occur within the original activator-inhibitor framework in the simplest setting in quantum dissipative systems. We note that the requirements of a coupled activator-inhibitor pair or the existence of homogeneous solution can be relaxed when we consider more general models [53][54][55][56][57][58][59] . In this study, we focus on the simplest case of a pair of symmetrically coupled quantum activator-inhibitor units and discuss quantum Turing instability in the original sense of Turing 23 and Gierer-Meinhardt 24 . Due to its simplicity, the model allows the direct numerical simulations of quantum dynamics and is the most amenable to experiment.
Correspondingly, in quantum system (3), when the diffusive coupling is absent ( D x = D p = 0 ), the state of each unit localizes around the stable fixed point at (0, 0) as shown in Fig. 2a. Thus, the two units obey the same distribution and the whole system is in the uniform state. However, when the diffusion constants are appropriately chosen, this uniform state is destabilized by the Turing instability and gives way to nonuniform states as demonstrated below. Figure 3 shows the Turing instability in the semiclassical regime observed by DNSs of QME (3). The same parameters as in Fig. 2 are assumed for both units. The two units are uncoupled ( D x = D p = 0 ) in Fig. 3a, 3c, 3e, while they are coupled with appropriate diffusion constants ( D x = 0.005, D p = 0.995 ) in Fig. 3b, 3d, 3f. To visualize the nonuniformity of the system state ρ , we introduce the two-mode Husimi Q distribution 72,73 Q x 1 , p 1 , x 2 , p 2 = 1 π 2 �α 1 , α 2 |ρ|α 1 , α 2 � with α j = x j + ip j (j = 1, 2) and use the marginal distributions Q(x 1 , x 2 ) = dp 1 dp 2 Q x 1 , p 1 , x 2 , p 2 and Q(p 1 , p 2 ) = dx 1 dx 2 Q x 1 , p 1 , x 2 , p 2 of the position (activator) variables x 1,2 and momentum (inhibitor) variables p 1,2 calculated from Q x 1 , p 1 , x 2 , p 2 .
In Fig. 3a, 3c without diffusive coupling, both Q(x 1 , x 2 ) and Q(p 1 , p 2 ) are symmetrically distributed around the origin. The variables of the two units are uncorrelated and statistically exhibit the same distribution. Thus, the state ρ of the whole system consisting of the two units is symmetric and uniform. In contrast, in Fig. 3b, 3d with diffusive coupling, Q(x 1 , x 2 ) is not symmetric and takes two extrema near the two classical fixed points (x 1 , x 2 ) = (A, −A) and (−A, A) , and similarly Q(p 1 , p 2 ) takes two extrema near (p 1 , p 2 ) = (B, −B) and (−B, B) . Thus, the two units tend to take the opposite states from each other and the state ρ of the whole system is nonuniform. It is noted that, because of quantum noise, the system state is mixed and the distributions have two symmetric peaks near both of the classical fixed points. Figures 3e and 3f show the marginal Wigner distributions W(x 1 , p 1 ) and W(x 2 , p 2 ) of units 1 and 2 for the cases without (e) and with (f) diffusive coupling. These Wigner functions are obtained from the marginal density operators ρ 1 = Tr 2 [ρ] and ρ 2 = Tr 1 [ρ] , where Tr j [·] represents the partial trace over system j in the semiclassical regime. Due to the symmetry of the two units, W(x 1 , p 1 ) and W(x 2 , p 2 ) are identical to each other. Additionally, the Wigner distributions in Fig. 3e without diffusive coupling are identical to that of a single unit shown in Fig. 2c. In Fig. 3e without diffusive coupling, the Wigner distributions have a single peak at the origin, whereas in Fig. 3f with diffusive coupling, the Wigner distributions have two symmetric peaks near the two stable fixed points (x 1 , p 1 , x 2 , p 2 ) = (±A, ±B, ∓A, ∓B) of deterministic classical system (4) (see "Methods").
The above results clearly indicate that Turing instability has indeed occurred and resulted in the formation of nonuniform stationary states in two diffusively coupled quantum activator-inhibitor units described by Eq. (3). In this regime, we can also perform direct numerical simulations of the corresponding SDE, which clearly visualize the nonuniformity caused by the Turing instability (see "Methods").
Weak quantum regime. Next, we show the results for the weak quantum regime. We set the parameters of QME (3) in a deeper quantum regime while keeping the deterministic system in the classical limit, Eq. (4), remain unchanged from the previous semiclassical case. See "Methods" for the characterization of the quantum www.nature.com/scientificreports/ regime. Figure 4 shows the Turing instability in this regime. The two units are uncoupled in Fig. 4a, 4c, 4e, while they are coupled with appropriate diffusion constants in Fig. 4b, 4d, 4f. As in the previous semiclassical case, when diffusive coupling is absent, the marginal Q distributions Q(x 1 , x 2 ) and Q(p 1 , p 2 ) of activator x and inhibitor p are symmetrically localized around the origin in Fig. 4a, 4c. When diffusive coupling is introduced, these joint distributions become nonsymmetric, indicating that the two units are anticorrelated and tend to take the opposite states from each other as shown in Fig. 4b, 4d. In this regime, due to the strong nonlinear damping, the two stable fixed points in the classical limit are closer to each other than in the semiclassical regime. Correspondingly, the nonuniformity of the joint distributions is less pronounced than in the semiclassical case due to the relatively strong effect of quantum noise.
Figures 4e and 4f show the marginal Wigner distributions W(x 1 , p 1 ) and W(x 2 , p 2 ) of units 1 and 2, which are identical to each other, before (e) and after (f) the Turing instability. Compared with the Wigner distribution in Fig. 4e before the Turing instability, the Wigner distribution in Fig. 4f after the instability is more elongated along the axis on which the two classical stable fixed points exist, although double symmetric peaks as in the semiclassical case are not observed due to the stronger effect of quantum noise.
Thus, although blurred by quantum noise, the system undergoes a transition from the uniform state to the nonuniform state with the introduction of diffusive coupling, namely, the Turing instability also occurs in the quantum regime considered here.

Strong quantum regime.
We also consider a strong quantum regime with a larger decay rate for nonlinear damping. Figure 5 shows the Turing instability in this regime. As the fluctuations are stronger than the two previous cases due to the effect of stronger quantum noise, only a slight nonuniformity can be observed. As shown later, the nonuniformity between the two units in this regime can be more clearly observed by using continuous measurement.
Phase diagram: nonuniformity and entanglement. We have seen that Turing instability occurs in a pair of diffusively coupled quantum activator-inhibitor units in the semiclassical, weak quantum, and strong quantum regimes. Here, we analyze the dependence of the system's behavior on the diffusion constants and the relationship between the Turing instability and quantum entanglement. We use the same parameter sets for the quantum activator-inhibitor units as in Figs. 3, 4, and 5 for the semiclassical, weak quantum, and strong quantum regimes, respectively. Figure 6 plots the (i) maximum eigenvalue max of the linearized equation of Eq. (4) in the classical limit (a, b),         Fig. 6a, 6b, the eigenvalue max of the uniform state is positive in the region below the dotted curve, where the diffusivity of the inhibitor D p is relatively large compared to that of the activator D x . Turing instability is expected to occur also in this region in the quantum system. The red dot ( D x = 0.005, D p = 0.995 ) represents the diffusion constants in the classical limit corresponding to Figs. 3, 4, and 5.
The RMSD plotted in Fig. 6c-6e shows that the nonuniformity is indeed caused by the Turing instability in the semiclassical, weak quantum, and strong quantum regimes and significantly correlated with the maximal eigenvalue max in the classical limit. There is a tendency that the nonuniformity is most strongly pronounced in the semiclassical regime (c), moderately in the weak quantum regime (d), and only weakly in the strong quantum regime (e), reflecting that the quantum noise is weaker and that the system state more clearly localizes around the two classical fixed points in this order (see Figs. 3, 4, and 5).
The negativity N shown in Fig. 6f-6h also increases with max , indicating that quantum entanglement between the two units also arises in the nonuniform state yielded by the Turing instability. Thus, the entanglement tends to be positively correlated with the nonuniformity between the two activator-inhibitor units and becomes stronger in the lower-right part where D x is small while D p is large in this parameter region. It is noted that a high-N region also arises when D p is close to zero while D x is relatively large, which is outside the Turingunstable region and simply shows that the two units are already entangled before the onset of Turing instability by the effects of two-mode squeezing and dissipative coupling.

Symmetry breaking via continuous measurement.
We have observed that Turing instability destabilizes the uniform state of the system of two units and gives rise to nonuniformity. The distributions in the nonuniform state are localized around the two classical fixed points as observed in Figs. 3, 4, and 5. This can be interpreted as a quantum-mechanically mixed state of the two classical situations where the system converges to either of the two stable fixed points. Thus, in contrast to the classical Turing instability in which only one of the two states is realized depending on the initial conditions, the symmetry of the coupled system is still preserved due to quantum noise even if the system state is nonuniform. Here, we show that further performing continuous measurement on the system can break this symmetry and reveal the true asymmetry of the system, which can www.nature.com/scientificreports/ be observed only in quantum systems. A similar measurement-induced spontaneous Z 2 symmetry breaking in a spin-chain system has been reported 74 .
We introduce continuous measurement on the linear damping (single-photon loss) bath coupled to each unit in QME (3). The stochastic master equations (SMEs) describing the system and the measurement results are then given by 75 where the first equation describes the stochastic evolution of the density operator ρ of the whole system under the effect of the measurement and the second equation describes the result Y j ( j = 1, 2 ) of the measurement on each unit. The term H[L]ρ = Lρ + ρL † − Tr [(L + L † )ρ]ρ represents the effect of measurement performed on the quadrature L + L † ; κ j and φ j (0 ≤ κ j ≤ 1, 0 ≤ φ j < 2π) represent the efficiency and quadrature angle of the measurement on the jth unit (j = 1, 2) , respectively; Y j is the output of the measurement result on the jth unit (j = 1, 2) ; and dW 1 and dW 2 represent independent Wiener processes satisfying �dW k (t)dW l (t)� = δ kl dt for k, l = 1, 2 . In contrast to QME, which gives averaged results over all possible measurement outcomes, this SME gives a single quantum trajectory of the system under the continuous measurement and can reveal the symmetry breaking of the system, which is preserved due to quantum noise in the steady state of QME. Figure 7 shows the behavior of the system under continuous measurement in the semiclassical regime. The parameters are the same as in Fig. 3b, 3d, 3f, namely, the uniform state of the system has been destabilized by the Turing instability. Considering that the nonuniformity is more pronounced in the position variable x than in the momentum variable p in Fig. 3d, we set φ j = 0 and perform the measurement on the quadrature x j = (a j + a † j )/2 ( j = 1, 2 ), which is conjugate to the momentum p j , of both units. We set the measurement efficiency as κ j = 0.25 ( j = 1, 2 ) for both units and the initial state of the whole system as the two-mode vacuum state. Figures 7a and 7b show the instantaneous marginal Wigner distributions W(x 1 , p 1 ) of ρ 1 and W(x 2 , p 2 ) of ρ 2 at time t = 50 sufficiently after the initial transient, obtained by a DNS of SME (5). In contrast to Fig. 3f, these Wigner distributions are not stationary and continue to fluctuate due to the continuous measurement. Each distribution is localized around either of the two stable fixed points of classical system (4) and tends to take the opposite state from the other one.
The anticorrelation between the states of the two units is evident in Fig. 7c-7f, where the time evolution of the average values of the position and momentum operators of both units, �x j � = Tr [((a j + a † j )/2)ρ] and �p j � = −iTr [((a j − a † j )/2)ρ] ( j = 1, 2 ), obtained from a single stochastic trajectory of quantum SME (5) are plotted. The two units randomly alternate between the two nonuniform states and tend to take opposite states from each other. This clearly indicates that the symmetry preserved by quantum noise is broken and that the asymmetry caused by the Turing instability in the classical sense is revealed by the extraction of information on the x variables of the two units via continuous measurement. Similarly, Fig. 8 shows the effect of continuous measurement in the weak quantum regime shown in Fig. 4. We observe qualitatively similar results to those for the semiclassical case in Fig. 7 in the quantum regime. Although the nonuniformity is less pronounced, the negativity is slightly larger on average, and the fluctuations are stronger due to the effect of the stronger quantum measurement noise. Notably, the negativity takes larger values than the case without performing measurement, indicating that the symmetry breaking due to the continuous measurement induces stronger entanglement in this regime.
Finally, we show in Fig. 9 the effect of continuous measurement in the strong quantum regime shown in Fig. 5. Although the fluctuations are stronger due to the effect of the stronger quantum measurement noise than the two previous cases, the nonuniformity between two single units, which was quite small in

Concluding remarks
We have theoretically demonstrated that Turing instability can occur in a quantum dissipative system. We showed that a degenerate parametric oscillator with nonlinear damping can be regarded as a quantum activator-inhibitor unit and that diffusive coupling between two such quantum activator-inhibitor units can give rise to Turing instability when the diffusivities of the activator and inhibitor variables are appropriately chosen. Due to the Turing instability, the system becomes nonuniform but still remains in a symmetrically mixed state by the effect of quantum noise. Further performing continuous quantum measurement breaks the symmetry and reveals the asymmetry between the two units. We suppose that the physical setup assumed in our model can, in principle, be implemented by using currently available experimental devices. The quantum activator-inhibitor unit is essentially a degenerate parametric oscillator with nonlinear damping 71 . The coupling terms via squeezing can be implemented by adjusting the single-mode squeezing parameter of the two quantum activator-inhibitor systems and introducing two-mode squeezing 76 . The dissipative coupling term could be realized by indirectly coupling the two oscillators through an additional cavity and adiabatically eliminating it 77 ; similar approaches have also been proposed for realizing dissipative couplings between ensembles of atoms 16 and optomechanical Stuart-Landau oscillators 14 . Another possible approach to the experimental realization of the proposed setups would be to use "membrane-in-themiddle" optomechanics 78 . Physical implementations of single-mode squeezing and nonlinear damping 79 , dissipative coupling 14 , and two-mode squeezing 80 have also been proposed. We expect that our numerical results for the Wigner distributions can be experimentally observed via quantum tomography 81 . The experimental implementation of the continuous quantum measurement has also been reported recently 82 .
In this study, we numerically analyzed a pair of quantum activator-inhibitor units that exhibits Turing instability in the classical, deterministic limit. For classical systems, analytical perturbative approaches have been applied to the classical master equation for predicting stochastic Turing patterns 41,83-85 . We may be able to employ similar perturbative approaches for the quantum master equation 12 and analyze the quantum Turing instability in more detail.
The quantum activator-inhibitor unit could also be implemented by using quantum spin systems, which is interesting because small quantum spin systems may help us cope with the exponential increase in the dimensions of the Hilbert space for large quantum networks 17 . Similar to previous studies that discussed the Kerr effects 15,86 and quantum jumps 87 in nonequilibrium pattern formation in quantum dissipative systems, clarifying the relationship between the Turing instability and strong quantum effects would be important. A more detailed systematic analysis on the relationship between Turing instability and entanglement is also a future study.
Although we analyzed only the minimal two-unit setup in this study, we may further consider Turing instability in larger networks of quantum activator-inhibitor units, similar to the Turing instability in networks of classical activator-inhibitor systems [45][46][47][48][49] . Compared to previous studies on quantum effects on nonlinear optical pattern formation 35,36 , which are not easy to analyze even numerically because calculations of all operator products are required 37 , the activator-inhibitor system proposed in this study can be extended to larger networks more easily. Thus, it may be used to reveal the novel emergence of self-organized patterns in quantum dissipative systems, similar to previous studies on the Kuramoto transition 12 , quantum chimera states 88 , and oscillation death 89 in globally connected quantum Stuart-Landau oscillator networks. Though we focused on a pair of coupled activator-inhibitor units in this study, we may also be able to further couple many units on a lattice or network of units and analyze the spatio-temporal pattern formation in fully quantum mechanical dissipative systems.
The quantum Turing instability may also find technical applications. For example, signal amplification near bifurcation points has been theoretically investigated in classical biological systems 90,91 and other classical 92 , nanoscale 93 , and quantum 94 nonlinear systems, and signal amplifiers using nonlinear bifurcation have been experimentally implemented 95 . Similarly, the Turing bifurcation in quantum dissipative systems may also offer new engineering applications for quantum signal amplification and quantum sensing.
As Turing instability is a paradigm of nonequilibrium self-organization in classical systems 96 , we believe that our results on the possibility of Turing instability in quantum dissipative systems also play an essentially important role in studying self-organization in quantum systems and will be relevant in the growing field of quantum technology.

Classical activator-inhibitor systems and Turing instability.
A classical activator-inhibitor system is generally described by where (˙) denotes the time derivative and x and p represent the activator and inhibitor variables, respectively. We assume that this system has a stable fixed point at (x, p) = (x,p) . Denoting small variations from (x,p) as δx = x −x and δp = p −p and linearizing Eq. (6), we obtain where we assume that the coefficients satisfy www.nature.com/scientificreports/ These are the conditions in which x is the activator and p is the inhibitor. These standard conditions can be eased in more general settings 55 , but we restrict our focus on the cases satisfying these conditions. We consider two diffusively coupled activator-inhibitor units with identical properties, described by where D x and D p represent the diffusion constants of the activator and inhibitor variables, respectively. This coupled system has a trivial fixed point (x 1 , p 1 , x 2 , p 2 ) = (x,p,x,p) , which corresponds to a uniform state of the whole system. In Turing instability, contrary to our intuition, this uniform state can be destabilized by the effect of diffusion when the parameters satisfy appropriate conditions. To see this, we linearize Eq. (9) as where δx j = x j −x and δp j = p j −p ( j = 1, 2 ) are small variations. The maximum eigenvalue of the Jacobian matrix in Eq. (10) is given by Therefore, when max > 0 , namely, when the uniform fixed point (x 1 , p 1 , x 2 , p 2 ) = (x,p,x,p) of the coupled system destabilizes.
In our model, the functions f and g are given by where γ 1 , γ 2 , η , and are parameters. The derivatives of f and g at this fixed point are given by With the parameter values used in the present study, the single system in Eq. (6) has a stable fixed point at (x, p) = (x,p) = (0, 0) , the conditions in Eq. (8) for the single system to be of the activator-inhibitor type are satisfied, and the condition for the Turing instability in Eq. (12) can be satisfied for a pair of diffusively coupled quantum activator-inhibitor units. As the Turing instability takes place, the trivial fixed point (0, 0, 0, 0) of the system is destabilized, and two new stable fixed points, which correspond to the nonuniform states of the whole system, arise via the supercritical pitchfork bifurcation, where (8) f x = ∂f /∂x| (x,p) > 0, f p = ∂f /∂p| (x,p) < 0, g x = ∂g/∂x| (x,p) > 0, g p = ∂g/∂p| (x,p) < 0. www.nature.com/scientificreports/ With the parameter values used in this study, the derivatives of f and g are f x = 0.5 , f p = −0.6 , g x = 0.6 , and g p = −0.7 . In Figs. 3 and 4, the maximum eigenvalue of the uniform fixed point is max ≈ 0.3724 > 0 ; hence, Turing instability has already occurred.
Quantum-classical correspondence via the Wigner distribution. We generally consider a quantum dissipative system with N modes, which is coupled with n reservoirs. We denote by a 1 , . . . , a N and a † 1 , . . . , a † N the annihilation and creation operators of the system, respectively. A general form of the QME describing this quantum dissipative system is given by where ρ is the density operator representing the system state, H is a system Hamiltonian, L j is a coupling operator between the system and jth reservoir (j = 1, . . . , n) , and D[L]ρ = LρL † − (ρL † L + L † Lρ)/2 is the Lindblad form 72,73 . By using the standard method of phase-space representation 72,73 , we can introduce the Wigner distribution W(α) ∈ R of ρ as and * indicates complex conjugate. QME (17) for the density operator ρ can be transformed into a partial differential equation for the Wigner distribution W(α) 72,73 , given by Here, the differential operator L p can be explicitly calculated from Eq. (17) by using the standard calculus 72,73 .
When the quantum effect is relatively weak, we may neglect the derivative terms higher than the second order in Eq. (19). Then, by introducing a real-valued representation of the phase-space variable, X = (x 1 , p 1 , . . . , x N , p N ) with α j = x j + ip j ( j = 1, . . . , N ), we can approximate Eq. (19) by the semiclassical FPE for W(X), Here, A(X) ∈ R 2N is the the drift vector, and D(X) ∈ R 2N×2N represents the diffusion matrix. The SDE corresponding to the above FPE is given by Here, A(X) is the same as in Eq. (20), the matrix G(X) ∈ R 2N represents the noise intensity satisfying G(X)G T (X) = D(X) with T representing the matrix transpose, and dW = (dw 1 , . . . , dw 2N ) ∈ R 2N represents a vector of independent Wiener processes satisfying �dw k (t)dw l (t)� = δ kl dt with k, l = 1, . . . , 2N . The deterministic trajectory in the classical limit is given by the deterministic term of the SDE, namely, Ẋ = A(X).
In the semiclassical regime where γ 2 is sufficiently small, the third-order derivative terms in Eq. (23) can be neglected 11,15,89 and the coefficients of the second-order derivative terms are positive. Therefore, Eq. (23) can be approximated by the FPE Using a real-valued representation, i.e., X = (x 1 , p 1 , x 2 , p 2 ) with α j = x j + ip j (j = 1, 2) , Eq. (25) can be rewritten as where Thus, the drift vector is given by A(X) = (A x 1 , A p 1 , A x 2 , A p 2 ) and the diffusion matrix D(X) is expressed as where we defined The SDE corresponding to FPE (26) is given by where G(X) satisfies G(X)G T (X) = D(X) and dW(t) = (dw 1 (t), dw 2 (t), dw 3 (t) , dw 4 (t)) T is a vector of independent Wiener processes satisfying �dw k (t)dw l (t)� = δ kl dt for k, l = 1, 2, 3, 4. When (30) dX(t) = A(X(t))dt + G(X(t))dW(t), www.nature.com/scientificreports/ Thus, the matrix G(X) can be chosen as G(X) = U (X) D ′ (X)U −1 (X) 89 , i.e., Direct numerical simulations of the quantum SDE. In addition to the QME, we also perform direct numerical simulations of semiclassical SDE (30) corresponding to FPE (26) to show the relationship of the distributions of the quantum states with the classical fixed points after the Turing instability. For example, supplementary Figures S1(a) and S1(b) show scatter plots of a stochastic trajectory of two diffusively coupled quantum activator-inhibitor units, and Fig. S1(c) shows the 2D plot of the Wigner distribution W(x 1,2 , p 1,2 ) in Fig. 3f. In Figs. S1(a) and S1(b), the states of units 1 and 2 stochastically go back and forth between the two stable fixed points due to quantum noise. These scatter plots agree with the Wigner distributions distributed around the two stable fixed points in Fig. S1(c).

Characterization of the quantum regime.
We characterize the degree of quantum effect as the nonlinear damping parameter γ 2 is varied by using the accuracy of the semiclassical approximation. The discrepancy between the semiclassical approximation and the original QME characterizes how deep the system is in the quantum regime. To keep the parameters of the corresponding classical systems unchanged, the linear damping parameter is chosen as γ 1 = γ ′ 1 + 2γ 2 , where γ ′ 1 is a constant, and the other parameters are fixed to the same values as those used previously. Figures 10a, 10b and 10c plot the average numbers of photons in both units and the nonuniformity �(x 1 − x 2 ) 2 � as functions of the nonlinear damping parameter γ 2 . Here, the average number of photons is calculated as an ensemble average �a † j a j � = Tr [a † j a j ρ] (j = 1, 2) of a † j a j obtained from the QME and as an average �α j α * j � α of α j α * j obtained from the semiclassical SDE, where the relation holds approximately in the semiclassical regime. The semiclassical results well approximate the results of the QME in the regime with small γ 2 , and the error due to the semiclassical approximation gradually increases with increasing γ 2 . Thus, when γ 2 = 0.1 (Figs. 2, 3, 6c, 6f, and 7), the semiclassical approximation is valid and the (34) (36) �α j α * j � α − 1/2 ≈ �a j a † j + a † j a j �/2 − 1/2 = �a † j a j � , weak quantum (f), and strong quantum regime (g). In (a-c), results obtained from the semiclassical SDE �α j α * j � α − 1/2 (red dots) and QME a † j a j (blue lines) ( j = 1, 2 ) are shown, where �α j α * j � α is calculated as a time average of α j (t)α * j (t) over a time interval of length 30000 after the initial transient. The parameters are � = −0.6, θ = π , η = 0.3 , D h = −0.99 , D c = 1 ( D x = 0.005 and D p = 0.995 ), and γ 1 = γ ′ 1 + 2γ 2 with γ ′ 1 = 0.2 . In (e-g), γ 2 = 0.1 (e), γ 2 = 0.5 (f), γ 2 = 3 (g). www.nature.com/scientificreports/ system is in the semiclassical regime, whereas when γ 2 = 0.5 (Figs. 4, 6d, 6g and 8) and γ 2 = 3 (Figs. 5, 6e, 6h and 9), the semiclassical approximation is no longer valid and the system is in the quantum regime. The degree of quantum effect can also be characterized by the purity as shown in Fig. 10d, where the purity increases with the increase of γ 2 . We also show in Fig. 10e-10g the elements of the density matrix of a single unit ρ 1 with respect to the number basis in the semiclassical (e), weak quantum (f), and strong quantum regime (g). We see that the energy level up to which the elements of the density matrix take non-zero value becomes lower and the discreteness of the energy spectrum becomes more prominent with the increase of γ 2 .
Negativity. We use the negativity N = ( ρ Ŵ 1 1 − 1)/2 to quantify the quantum entanglement of the two units, where ρ Ŵ 1 represents the partial transpose of the density operator ρ of the two-mode system with units 1 and 2 with respect to unit 1 and �X� 1 = Tr|X| = Tr √ X † X 97,98 . A non-zero negativity indicates that the two units are entangled. Note that the negativity N ′ = ( ρ Ŵ 2 1 − 1)/2 calculated with respect to unit 2 is equal to the negativity N calculated with respect to the unit 1.

Data availability
All data generated or analysed during this study are included in this published article and its supplementary information files.